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1. Motivation and objectives 

Turbulent reacting flow has been an important problem and has attracted much 
attention from researchers in a variety of science and engineering disciplines. Despite 
intense research activity, however, much remains to be done in this field (O’Brien 
1980; Pope 1985, 1990). One of the key issues in engineering application is to 
employ the existing models and techniques to develop a relatively simple numerical 
scheme for simulating complicated reacting flow systems. 

Several approaches have been introduced to overcome the closure problems en- 
countered in turbulent reacting flow simulations. Among them, the probability 
density function (PDF) method provides a closed form representation for the chem- 
ical source terms (O’Brien 1980; Pope 1985); for this reason it has become a pre- 
ferred choice. However, the scalar PDF does not contain information concerning 
the transporting velocity field and the interactions between the scalar and velocity 
fields that has to be provided by supplementary turbulent transport and mixing 
models in reacting flow simulations. 

The recent development of the dynamic subgrid-scale (SGS) model (Germano 
1992; Germano et al 1991; Ghosal et al 1992) has provided a consistent method for 
generating localized turbulent mixing models and has opened up great possibilities 
for applying the large eddy simulation (LES) technique to real world problems. 
Given the fact that the direct numerical simulation (DNS) can not solve for engi- 
neering flow problems in the foreseeable future (Reynolds 1989), the LES is certainly 
an attractive alternative. It seems only natural to bring this new development in 
SGS modeling to bear on the reacting flows. 

The major stumbling block for introducing LES to reacting flow problems has 
been the proper modeling of the reaction source terms. Various models have been 
proposed, but none of them has a wide range of applicability. For example, some 
of the models in combustion have been based on the flamelet assumption (Kerstein 
et al 1988; Trouve & Poinsot 1992), which is only valid for relatively fast reactions. 
Some other models have neglected the effects of chemical reactions on the turbulent 
mixing time scale (Valino Gao 1992), which is certainly not valid for fast and 
non-isothermal reactions (Vervisch 1993). 

The PDF method can be usefully employed to deal with the modeling of the 
reaction source terms. In order to fit into the framework of LES, a new PDF, the 
large eddy PDF (LEPDF), is introduced. This PDF provides an accurate represen- 
tation for the filtered chemical source terms and can be readily calculated in the 
simulations. The details of this scheme are described below. 
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2. Accomplishments 

Large eddy PDF 

The large eddy fields, which are explicitly simulated in the LES, can be obtained 
by filtering the true fields with certain filters G (Germano 1992; Rogallo k Moin 
1984); namely, 

A(x, t)= f A(x', t)G(x' - x)dx'. 

J — OO 

Among the commonly used filters, we are particularly interested in those that are 
localized in physical space, such as the local volume average (Schumann 1975) and 
the Gaussian filters (Leonard 1974), since they describe local averaged effects. For 
reasons that will become clear later, we choose only these positive definite filters. 

By applying a filter of size A, which is generally the mesh size in LES, the Navier- 
Stokes equation can be written as 


du, _ dui 
dt + U *dxj 


— — Vp — 



( 1 ) 


Here r tJ = uTuj - u,u } is the SGS stress and is normally modeled by the eddy- 
viscosity model originally proposed by Smagorinsky (Smagorinsky 1963): 


T ij - \r kk 6 tj = — 2CA 2 |S|3>,-j, 


( 2 ) 


where S tJ is the strain rate tensor 

5 - = 1(55 i + «5i) 

^ 2 ( dxj + dxi ] 


and \S\ = ^2S t jS l} . 

The same type of filter with larger size A > A can be applied to the same 
equation. The resulting SGS stress can be represented by 

Tij = UiUj - UiUj. 

If the filters are well behaved ones such as the Gaussian filters, we will have the 
following convolution relation 

G(x — x',l) = G(h)* G(h) = J G(x — xi,li)G(xi — x' ,h)dx\. (3) 

Take the Gaussian filter as an example; it can be shown that l 2 = l\ + 1%. It is 
obvious from (3) that A - % where “A” is the filtering performed by the “gap 
filter” between A and A. The Germano identity (Germano 1992) 

T,j - Tij = Lij = UiUj - uitij, 


( 4 ) 
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where the right hand side can be explicitly evaluated from the LES field, can then 
be used to calculate the local coefficient C for the Smagorinsky model (Germano et 
al 1991; Moin 1991; Lilly 1992). 

A similar method can be applied to the turbulent scalar field which is governed 
by 

dib dib 0 

(5) 

Its large eddy counter-part can be written as 


dib dib o — 9 

+ u jfc~ = D* 2 * ~ ^-(w>V’-«>V’) + w(x,<). (6) 

Similar eddy- viscosity type models can be developed for the scalar field (Moin et al 
1991), namely 

_ — dib 

ujtb - ujib = ( 7 ) 

where D a can again be obtained through the dynamic procedure described above. 

The major problem is to find a suitable SGS model for the reaction source term 
w(^). Noticing that 

<^( x >*)=/ u(xl>(x',t))G(x' -x)dx' = f[d<j>u>(4>) f p(<t>;x' ,t)G(x' - x)dx'\, 

J — OO J J — OO 

( 8 ) 

where 

p(<f>;x,t) = 6(<f> - xl>(x,t)) 

is the fine-grain function, it is clear that once 


Pl(4>; x, t) = f p{<f>\ x', t)G(x' - x)dx' 

J — OO 


( 9 ) 


is known, the filtered equation for the scalar field %b I s closed. For positive- definite 
filters, it can be shown that Pi has all the properties a PDF should have. We define 
it as the large eddy PDF. It is a generalization of the existing PDF concept. 

In fact, for a homogeneous field, where the ensemble average can be replaced by 
the space average, the traditional PDF (O’Brien 1980; Pope 1985) can be recovered 
from our definition if a volume average filter is applied to a large enough space 
volume. For more general cases, Pi can be regarded as a weighted average of all 
contributions from the neighboring field points, depending on how far they are from 
the observation point. 

The governing equation for Pi can be derived following the standard procedure 
(O’Brien 1980; Pope 1985). It can be shown that 
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Applying the filter on both hand sides of (10), we get 

^ + u • = -^[dW^ + u^Pl] - ^-(tiT p-Ujp). (11) 

The last term on the RHS of (11) can again be modeled by the eddy viscosity 
assumption. Since p is controlled by the underlying process ip, the D, obtained for 
ip is equally applicable for p (Jiang & O’Brien 1991). A brief proof follows. 
Suppose the eddy- viscosity for Pi is D p . Note that 

m=f <ppu<p)d<p. 

It can be shown, by multiplying <f> on both sides of equation (11) and integrating 
over the composition space after inserting 

__ _ _ n dP L 

UjP- u jP= - D p-fa~ 


that 




Comparing this equation to equations (6) and (7), it is obvious that D p 
It can be shown that 


= D a 


(pV*ij>) = E L {V 2 ii>\<f>}P L (<f>). 

If A and B are composition space representation of the fields a(x) and fr(x), then 
Ei is defined by 

E l {A\B}P l (B) = J APi(A,B)dA 

and 

Pl(A,B) = I 6(A- a(x'))6(B - t(x'))G(x' - *)dx'. 

Ei{A\B } can be interpreted as the average of A, weighted by G(x' — x), over 
the spots where b(x') = B . The modeling of this term is generally difficult, as 
we have experienced in the traditional PDF approach. However, the fundamental 
physics expressed by this term remains the same - it represents the enhancement 
of diffusion by turbulent fluctuations. Therefore, models analogous to those used 
in the traditional PDF can be used to close this term. For example, if the LMSE 
model (O’Brien 1980) were used, we would have 

E L {V 2 ip\<p} ^-—(<p-4,) + dv 2 4>, 
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where I/77, oc ( D + D s ). After substituting this relation into (11), it can be shown 
that (11) recovers (6). Therefore, this relaxation model is adequate if one is only 
interested in the large eddy field. 

Equation (11), supplemented by the above closure models, can be solved using 
the well established Monte-Carlo techniques for the scalar PDF (Pope 1981, 1985). 
A few advantages of this simulation scheme can be immediately identified. First, a 
general model for the chemical source term is developed based on the PDF. It, in 
principle, applies to all kinds of reactions and does not require any artificial pre- 
assumptions. Secondly, the interactions between chemical reaction and turbulent 
mixing are embedded in the current model - the effects of reaction on turbulent mix- 
ing are reflected in D s through the dynamic procedure. Thirdly, the LES provides 
the flexibility of resolving a certain range of large scale structures while modeling 
the contributions from the rest of the scales. It can be applied to a wide range of 
problems with different degrees of resolution depending on the available resources. 
Moreover, it may be worth mentioning that the Monte-Carlo simulation for scalar 
PDF is a well established technique while LES has shown great potential with the 
introduction of the dynamic SGS model. The outlook for applying the proposed 
scheme seems to be promising. We intend to apply the scheme described above 
to reactive flow problems in a number of turbulent flows, such as in homogeneous 
flows, for which DNS results are available. 

A consistency condition for SGS model 

There are, however, certain problems that should be further investigated. Cor- 
rectly modeling turbulent mixing in the PDF formulation has been a long-standing 
problem which surely deserves more attention (Gao 1993). The consistent determi- 
nation of the eddy viscosity in the LES is another such problem. 

Several different consistency conditions should be considered in calculating the 
eddy viscosity. First, the eddy viscosity is localized; namely, it is a function of both 
space and time. It can not be taken freely out of the space filter; therefore, it cannot 
be evaluated by simple algebraic procedures. A variational approach (Ghosal et al 
1992), which leads to an integral equation for the coefficient C, and a differential 
equation approach, which will be described below, have been developed to overcome 
this problem. Another problem is that there is no reason to believe that the eddy 
viscosity should be the same on both the cut-off (A) and the test (A) filter levels. 
Since the equation for C constitutes an over- determined problem (six equations for 
two C in the case of velocity, if different model coefficients axe assumed for different 
filter levels), the least mean square estimate technique can be used to evaluate these 
different model coefficients at different levels (Moin 1991; Lilly 1992). To use the 
variational approach, two integral equations will have to be solved. It could be 
numerically very intense to do so. 

In fact, if A and A are in different scale ranges, there is no reason even to believe 

that Tij ~ Tij(S t j) and r t j ~ T ij(Sij) will have the same functional form. It is, 
therefore, helpful to imagine that the cut-off and the test levels are brought closer 
and closer to each other. In this limiting process, we can be assured that T{j and 
Tij will have the same functional form and that the coefficients C at different levels 
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will approach each other. 

A formal procedure can be performed in this limiting case. For simplicity, we will 
omit the bar in u in the following derivations. 

We assume that the gap between the cut-off and the test filters Eire represented 
by e. As the gap between two filters narrows (e — ► 0), G(x' - x, «) becomes more 
“concentrated”. In fact, 

lim G(x' — x, e) = 6(x! — x). 

t — ►O 

Let 

9 = J(xi - x) 2 G(xi - x,e)dx!, (12) 

it is obvious that lim £ _^o 0 = 0. 

For simplicity, we will only consider symmetric filters, i.e. G = G( |xj — x|). 
However, this method is not restricted to symmetric filters. 

Under symmetric conditions, we have 

u(x) — J uG(x i - x, e)dx i = u(x) + ^ V 2 u + O(0 2 ). (13) 

Obviously, 

Lij = 6Vui * Vuj + 0(0 2 ). (14) 

Similarly, for the Smagorinsky type of SGS model, we have 

T„ - t,j - + !|!v 2 S„ + 

+ A^(C|S|S„) + |^|t^] + 0(^), (15) 

where we have used the relation A 2 = A 2 + a0 2 . For Gaussian filters, a = 1. 

Substituting both (14) and (15) into the Germano identity (4) and taking the 
limit e — ► 0, we have 

-2CA 2 ( 0 ^Sii + — Sjj) + A 2 V 2 (C|S|5„) 

-Wm-V., (16, 

This equation, when properly contracted, can be used to solve for the local model 
coefficient, C . Clearly, this procedure can be equally useful if different functional 
forms for Sij and are chosen. 
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